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BUOYANCY  DRIVEN  FLOW  AS  THE  FORCING  FUNCTION  OF  SMOKE  TRANSPORT  MODELS 


Walter  W.  Jones 
Xavier  Bodart 

Abstract 

Flow  at  vents  is  the  major  driving  force  in  smoke  transport  models.  The 
precision  with  which  we  can  calculate  these  flows  determines  to  a great  extent 
how  accurately  we  can  model  buoyant  flow  and  the  inherent  speed  of  the  models. 
This  report  describes  some  of  the  problems  encountered  in  calculating  these 
flows,  and  gives  a general  algorithm  for  their  calculation. 

Key  words:  vent  flow,  smoke  transport  model,  fire  modeling,  differential 

equation. 


1 . INTRODUCTION 


Over  the  past  few  years  modeling  fire  growth  and  smoke  transport  has 
become  an  important  aspect  of  fire  research.  As  computers  have  become  cheaper 
and  faster,  the  ability  to  handle  the  equations  associated  with  such  phenomena 
has  improved  to  the  point  where  numerical  experiments  can  replace  some 
physical  experiments.  This  is  useful  in  avoiding  the  expense  and  effort  of 
actual  fire  tests.  To  this  end  there  are  two  aspects  of  modeling  which  are 
important.  The  first  is  that  the  physical  algorithm  must  be  correct  and  the 
second  is  that  the  numerical  routines  which  are  utilized  must  be  fast  and  able 
to  handle  the  wide  range  of  values  which  occur  in  natural  phenomena.  Typical 
time  scales  range  from  microseconds  for  chemical  kinetics  to  tens  of  seconds 
for  heat  conduction. 
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The  predictive  equations  which  are  solved  arise  from  the  conservation  of 
mass,  momentum  and  energy.  When  integrated  over  a finite  sized  control 
volume,  they  are  first  order,  nonlinear,  ordinary  differential  equations,  of 
the  form 

^ = f(x),  (1) 

where  x represents  a vector  of  unknown  dependent  variables.  Much  of  the 
modeling  has  assumed  a simplification  which  is  to  solve  the  pseudo  steady- 
state  as  the  transient  term  vanishes 

f -£(x)  >0. 

In  either  case  a solution  is  found  by  varying  x in  the  phase  space  about  some 

initial  value  x and  looking  for  a minimum  in  the  value  of  the  function  f(x) 
o 

or  actually  - f(x)}.  The  thesis  of  this  paper  is  the  difference  in  the 
physical  meaning  of  these  two  types  of  equations  and  to  demonstrate  that 
solving  the  original  ODE's  is  superior  to  solving  their  algebraic  counterpart. 
The  underlying  point  is  to  demonstrate  the  rationale  that  the  apparently  more 
complex  form  of  the  equations  is  actually  easier  to  solve,  and  in  addition  to 
show  that  the  algebraic  form  may  not  always  yield  the  correct  solution. 

2.  PROBLEM  TO  BE  SOLVED 

The  difficulty  which  arises  in  the  solution  of  eqn.  (1)  is  the  bumpiness 
in  the  multidimensional  phase  space  allowed  by  these  equations  (one  for  each 
zone  or  control  volume).  The  term  "phase  space"  is  used  in  this  context  to 
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mean  the  manifold  of  pressure,  temperature,  etc.,  in  which  the  solution  of  the 

relevant  set  of  equations  lie,  and  the  possible  physical  variations  which  can 

occur.  Even  for  a single  equation,  e.g.  the  pressure  equation,  the  topology 

is  not  simple.  There  is  an  assumed  simplification  of  considering  only  relaxed 
dx 

states  — = 0 which  can  actually  make  achieving  the  final  state  more  diffi- 
dt 

cult.  No  direction  is  provided  in  finding  a minimum  and  it  gives  no  sense  of 
whether  a specific  bump  in  the  manifold  is  an  absolute  minimum.  By  retaining 
the  transient  term  some  relief  is  obtained  for  this  situation.  The  most 
important  improvement  in  finding  a solution  is  that  search  direction  in  phase 
space  is  given  by  the  transient  term.  Thus  the  basic  search  algorithm  is 
considerably  simplified.  The  second  bit  of  help  comes  from  the  observation 
that  the  derivative  terms  do  not  vanish  at  a pseudo  minimum  but  only  at 
absolute  minimum.  Further,  the  point  for  which  the  transient  term  vanishes 
may  not  be  the  correct  solution.  The  correct  solution  actually  minimizes  the 
difference  between  the  left  and  right  side  of  eqn.  (1). 

Several  examples  should  serve  to  illustrate  the  problem,  a possible 
solution  and  the  means  by  which  this  translates  into  a technique  for  a general 
solution  of  the  equations. 

The  actual  equations  which  we  solve  are  somewhat  more  complex  than 
eqn.  (1);  however,  it  is  useful  to  use  a simplified  example  to  test  the 
techniques.  Then  one  can  test  the  heuristic  without  getting  bogged  down  in 
technical  details  of  numerical  programming.  The  set  of  predictive  equations 
which  we  use  for  each  compartment  in  a typical  zone  model  [1],  are 


dt 


(3-l)V 


s 


(2) 
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dT 


dt 


1 / +— S—s 


® “ (6-1) 


dt  3 \PV^j(^  A 


dV 
u 

dt 


■(h) 


V 

c m T + E - — ^ s 
p u u u ^ 


where 


s = c m T + c m.T.  + E + E. 
puu  p Z Z u Z 


and 


3 = c /R  = y/(y-1) 
P 


(3) 

(A) 

(5) 


where  T is  the  temperature  and  E^  and  E^^  are  the  energy  release  rates  into  the 
upper  and  lower  zones,  respectively.  The  primary  source  term  is  the  s which 
depends  mainly  on  fluid  transport,  that  is  enthalpy  flux  driven  by  density  and 
pressure  gradients.  It  also  varies  the  most  rapidly.  Actually,  the  micro- 
scopic chemical  kinetics  have  a shorter  time  constant  but  the  spatially 
averaged  zone  properties  (T,  p)  do  not  vary  on  such  a time  scale  and  thus  it 
is  reasonable  to  ignore  these  very  short  times.  In  addition  to  improving  the 
solution  technique  itself,  we  have  developed  a better  method  of  finding  the 
dominant  portion  of  the  source  term,  m,  the  enthalpy  flux  from  mass  flow.  Due 
to  its  importance,  we  will  discuss  it  in  the  next  section  before  giving 
examples  of  the  solution  of  this  type  of  problem. 


The  conservation  equations  which  we  solve  are  for  energy  and  mass. 
Momentum  is  not  defined  within  a control  volume  but  only  at  the  boundaries 
where  zones  are  connected.  So  we  do  not  solve  the  momentum  equation  directly 
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but  rather  use  the  integral  form  known  as  Bernoulli’s  equation  which  yields 
mass  flux  as  a function  of  the  pressure  and  density  differences. 


3 . SOURCE  TERM 

Equation  (2)  is  the  one  on  which  we  shall  focus  as  it  is  this  one  that  is 
most  often  "simplified"  in  the  as3nnptotic  sense.  In  order  to  explain  some  of 
the  difficulty  encountered  in  solving  this  equation,  we  digress  for  a moment 
to  discuss  the  primary  source  term,  namely  fluid  transport  through  vents. 

This  fluid  flow  phenomenon  connects  the  control  volumes  and  is  dominant 
because  this  term  fluctuates  most  rapidly  of  all  the  source  terms  in  response 
to  changes  in  the  environment.  One  of  the  inq)rovements  which  we  have  incorpo- 
rated into  our  current  models  is  a means  of  calculating  these  flow  fields  with 
the  correct  number  of  neutral  planes  (up  to  three)  and  without  discontinuties 
in  the  function.  This  latter  feature  implies  no  discontinuities  in  the  first 
order  derivatives  for  the  ODE's. 

Typical  types  of  flows  which  can  occur  in  fires  are  illustrated  in 
fig.  1.  The  notation  for  the  flow  is 

SS  = upper  layer  to  upper  layer 
SA  = upper  layer  to  lower  layer 
AS  = lower  layer  to  upper  layer 
AA  = lower  layer  to  lower  layer 
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This  notation  was  originated  earlier  [2]  for  the  single  neutral  plane  case  but 
is  useful  for  a physical  description  of  the  general  problem.  One  generally 
uses  the  Bernoulli  equation  to  calculate  the  flow  velocities  between  two 
compartments  which  are  connected  by  an  opening.  Indeed  this  is  the  solution 
for  the  momentum  equation  which  allows  us  to  exclude  it  specifically  when 
solving  the  conservation  equations  in  general. 

The  general  form  is 


• 

m.  = 

lO 

C 

• S • /2p  (kg/sec). 

where  m 

= 

mass  flow  rate 

C 

= 

orifice  coefficient 

S 

opening  area  (m^) 

P 

= 

gas  density  on  side  "i" 

P . 
1 

= 

pressure  on  side  "i" 

P 

o 

= 

pressure  on  side  "o". 

The  implication  of  using  this  equation  is  that  the  pressure  at  a stagnation 
point  is  used.  That  is,  the  flow  velocity  vanishes  where  the  pressure  Is 
measured.  As  the  pressure  always  appears  as  a difference,  in  principle  it 
does  not  matter  whether  absolute  pressures  or  pressure  defects  are  used.  More 
accuracy  is  obtained  by  using  pressure  defect,  however.  This  avoids,  to  some 
extent,  the  problems  of  the  small  difference  of  large  numbers. 

The  pressure  is  always  calculated  with  respect  to  the  base  of  a compart- 
ment. With  this  in  mind  we  can  express  the  pressure  on  the  other  side  of  a 
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partition  as  a function  of  the  variable  (y),  the  height  above  the  base: 
implicit  in  our  use  of  equation  (6)  is  that  the  opening  is  rectangular,  so 
that  the  area  integral  of  the  flow  term  will  allow  us  to  remove  the  width  from 
the  integral. 

That  is 

Z2 

flow  = / /,  . , . pV  dzdb  width  / pVdz.  (7) 

•'width  ■’height  Z' 

Thus  the  total  flow  becomes 

m_  = C W Z S,  /2  • p.  , |p.(z)  - P (z)|  . 

i->-o  ^ k l»k  I 1 o ' 

The  pressure  term  will  be  reversed  if  the  flow  is  o ->■  i.  Thus  we  have  the 

integral  over  the  area  shown  in  fig.  2.  Terms  are  as  before  except  that 

p.  , is  the  average  inlet  mass  density  within  area  "k".  The  simplest  way  to 
1 , K 

define  the  limits  of  integration  is  with  neutral  planes,  that  is  the  height  at 

which  flow  reversal  occurs,  P.  (z)  = P (z). 

1 o 

Each  side  of  an  opening  is  assumed  to  consist  of  two  homogeneous  gas 
layers  (zones)  of  uniform  density  and  temperature.  There  is  an  apparent 
inconsistency  in  that  the  equation  of  state  dictates  P = pRT,  and  we  assume  P 
varies  but  p and  T remain  constant,  at  least  within  a zone.  This  pressure 
fluctuation  is  so  small  compared  to  the  magnitude  of  the  base  pressure  that 
ignoring  it  for  all  calculations  except  the  flow  field  is  reasonable.  For 
simplicity  of  notation  we  will  use  a slightly  different  means  of  Identifying 
the  zones.  The  correspondence  is 
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upper  "i"  zone  = 1 
lower  "i"  zone  = 2 
upper  "o"  zone  = 3 
lower  "o”  zone  = 4 

The  former  terra  using  "SA"  is  physically  more  understandable  and  the  results 
can  be  put  into  these  terms,  but  the  derivations  are  more  compact  with  the 
numerical  indicies.  Figure  3 shows  a schematic  of  the  notation.  The  corres- 
ponding densities  are  , P^  p^.  and  are  the  height  of 

a sill,  soffit,  hot/cold  interface  in  the  "i"  compartment  and  hot/cold  inter- 
face in  the  "o"  compartment.  With  the  base  (reference)  pressures  P^(0)  and 
Pq(0),  the  internal  pressure  on  each  side  is  given  by 

P^  (z)  = P^(0)  - min  (y,  Z^)  • p^g  - max  (z-Z^,  0)  • p^g 

P (z)  = P (0)  - min  (y,  Z ) • p, g - max  (z-Z  , 0)  • p g 
o o o 4 o J 

The  function  F (z)  = P. (z)  - P (z)  can  be  considered  a family  of  functions  of 

i e 

one  variable,  z.  In  all  cases  the  pressure  appears  only  in  this  form.  In 
principle,  this  family  of  curves  can  lead  to  an  inordinate  number  of  possible 
flow  fields.  By  Imposing  the  restrictions  found  in  fire  scenarios,  however, 
we  end  up  with  only  a few  possiblities . For  stratification  to  occur,  the 
following  restrictions  can  be  imposed; 


We  can  also  require 


without  loss  of 
reverse  the  "i" 
These  cases  and 


generality,  since  for  the  reverse,  it  is  only  necessary  to 
and  "o"  compartments.  We  are  left  with  five  different  cases, 
their  restrictions  are  shown  in  Table  I. 


TABLE  I 


Class  Restrictions  Max  //  of  neutral  planes  Figure  4 


1 

^2 

Q. 

z . 
1 

< z 

— 0 

1 

a 

II 

^2 

< 

a. 

z . 

1 

< z 

— 0 

2 

b 

III 

^3 

< 

^2 

< 

P, , Z . > z 

4 1 0 

3 

c 

IV 

^2 

< 

0K 

Cl 

z, 

1 

< z 

0 

2 

d 

V 

^2 

< 

P3. 

z . 
1 

< Z 

0 

1 

e 

Classification  of  types  of  flow  which  occur  in  a vent  based  on  the  relative 
densities  and  interface  heights.  If  there  exist  soffits  or  sills,  then  the 
number  of  neutral  planes  within  a vent  can  be  less  than  the  number  indicated. 


If  there  were  no  soffits  or  sills  to  consider,  then  the  calculation  would  be 
fairly  straightforward.  However,  the  possibility  of  soffit/sill  combinations 
requires  many  numerical  tests  in  the  calculation.  Class  I is  the  basis  of  the 
analysis  of  classes  11-V.  It  can  have  44  different  flow  combinations,  depend- 
ing on  the  relative  position  of  H^,  B^,  and  Z^.  It  may  contain  at  most  a 
single  neutral  plane  (flow  reversal).  Twenty  four  of  these  combinations  are 
without  a neutral  plane  and  twenty  with  a neutral  plane.  Figure  5 shows  the 
effect  of  the  "o"  compartment  interface  height  on  the  flow  field  from  one  of 
the  diagrams  shown  in  Figure  4. 
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Table  II  shows  the  criteria  hsed  for  solution  in  classes  II-V.  The 
interval  H^]  can  subsequently  be  partitioned  to  contain  at  most  a single 

neutral  plane  so  that  the  logic  used  for  class  I can  be  utilized.  The  only 
caveat  is  to  be  sure  that  the  equalities  and  inequalities  are  treated 
correctly. 


4.  A SAMPLE  PROBLEM 


We  limit  our  discussion  to  eqn.  (2)  with  the  appropriate  source  terms. 
It  is  most  often  this  equation  to  which  the  asymptotic  approximation  is 
applied,  namely 


dt 


-►  0. 


(8) 


This  approximation  is  appropriate  for  a steady  state  situation,  for  example, 
where  a fire  is  fully  developed  and  the  flow  fields  have  been  established. 

For  a transient  problem  such  as  a developing  fire  or  the  case  when  a window  Is 
broken,  such  is  not  appropriate.  Most  models  of  fire  imply  this  approximation 
by  using  large  initial  time  steps,  perhaps  of  one  or  two  seconds.  However, 
when  faced  with  a true  transient  phenomenon,  all  sub-divide  the  time  step, 
generally  to  a value  below  the  pressure  relaxation  time. 


Furthermore,  a fundamental  problem  arises  in  that  the  solution  for 
eqn.  (2)  may  not  be  identical  to  that  of  eqn.  (8).  The  requirement  Is  that 
the  difference  between  the  left-  and  right-hand  sides  (LHS  and  RHS)  of  the 
conservation  equations  should  be  minimal.  This  is  truly  equivalent  to  the 
statement  of  eqn.  (8)  only  in  a steady  state  regime.  So  this  approximation 
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TABLE  II 


F(z)  = P^(z)  - Pg(z) 


Condition  Max.  no.  of  neutral  planes 

in  [Bj,  H^] 


II.  j<  B^ 


0,  1 


Zi  > Hf 


Bf  < Z.  < 

F(Zi)  • max  (F(Bf),  F(Hf))20 
F(Z^)  • max  (F(Bj),  F(H^)  < 0 

F(Zj^)  • min  (F(Bf),  F(H^))  2 0 
F(Z^)  • min  (F(Bf),  F(Hf)  ^0 

III.  2 Hf 


Zo  <.  Bf,  Z^  2 Hf 


2 Bf, 


Zi 


< Bf 


0,  1 

0 

1 

2 

0,  1 

0 
0 


Zo  < Bf* 


Bf  < Zf  < Hf 


F(Zf)  • max  (F(Bf),  F(Hf))  2.0  0 

F(Z^)  . max  (F(Bf),  F(Hf))  > 0 

F(Zf)  • min  (F(Bf),  F(Hf))  20  1 

F(Z^)  • min  (F(Bf),  F(Hf))  2 2 


B 


f 


< Zo  < Hf 

Zi  > Hf 

F(Zq)  • min  (F(Bf),  F(H.))2  0 0 

F(Z^)  • min  (F(Bf),  F(Hf))  < 0 

F(Zq)  max  (F(Bf).  F(H.))  1 

F(Z  ) • max  (F(B.),  FtHf))  2 0 2 

Z^  < H^  [Bj  < Z < Z,  < Hj^] 

min  (F,(Bf1,  F(Z.))  • max  (F(Hf),  F(Z  ))  >00 

min  Fj(B^),  F(Z^7)  * max  (F(Hf),  F(zJ)  >_  0 1,2  or  3 


IV.  Z . 2 Hf  0,1 

Zf  2 Bf  0,1 

Bf  < Z.  < Hf 

FCZf)  • max  (F(Bf),  F(Hf))2  0 0 

F(Zf)  • max  (F(Bf),  F(Hf))  < 0 

F(Zf)  * min  (F(Bf),  F(Hf))  > 0 1 

F(Z.)  • min  (F(Bf),  F(Hf))  2 0 2 


V.  F is  monotonic  over  [B^,  Hf] 


0,1 
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should  only  be  made  if  there  is  a real  gain,  such  as  reduced  computing 

I 

requirements  or  equations  which  are  simpler  to  manipulate.  An  example  should 
Illustrate  some  of  the  difficulties. 

The  context  of  a problem  which  we  will  examine  is  a four  compartment 
calculation  using  the  FAST  [1]  model  to  find  a solution.  The  point  of 
Interest  is  to  examine  the  solution  of  eqn.  (2)  in  the  second  compartment. 
There  are  three  ways  to  examine  the  equation;  first,  a pseudo  analytic  tech- 
nique is  to  start  with  the  converged  solution  of  temperature  and  pressure  in 
all  compartments  and  then  form  a perturbation  expansion  of  the  RHS  of 
eqn.  (2),  in  terms  of  pressure  and  temperature  for  compartment  //2  only.  The 
result  is  shown  in  fig.  (6)  as  a surface  plot  of  the  value  of  the  PHS  of  eqn. 
(2)  for  a variation  in  pressure  (left  to  right)  and  temperature  (into  the 
picture).  The  total  calculation  lasts  60  seconds  and  fig.  (6a,  6b,  6c,  and 
6d)  are  done  at  0,  20,  40,  and  60  seconds,  respectively.  The  variation  in 
pressure  was  temperature  T^  100  K.  If  this  equation  were 

solved  by  Itself,  then  the  solution  would  come  fairly  easily  and  eqn.  (8) 
might  be  appropriate;  second,  we  can  capture  the  results  from  the  ODE  solver 
which  is  used  in  FAST.  This  was  done  at  50  ± 1 seconds  and  the  results  are 
shown  in  fig.  (7).  The  axes  are  the  same  as  in  fig.  (6).  The  latter  figure 
shows  the  somewhat  more  complex  interaction  of  all  sixteen  equations,  since  we 
capture  the  intermediate  results  with  variations  in  all  parameters.  From  the 
figure  it  can  be  seen  that  it  may  not  be  possible  to  get  from  a pseudo  minimum 
on  the  left  hand  side  to  the  correct  answer  on  the  right  hand  side  using 
eqn.  (8).  Equation  (2)  helps  by  providing  a direction  and  distance  for  subse- 
quent solutions,  namely  the  dP/dt  term,  the  term  which  has  been  dropped  In  thi- 
asymptotic  approximation.  Figure  (8)  shows  a composite  of  the  result  of 
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applying  the  second  technique  to  the  entire  time  history  of  the  pressure 
equation. 


Finally  we  can  look  at  a simple  piece  of  this  problem  analytically.  We 
can  extend  earlier  w/ork  on  pressure  relaxation.  Attention  is  directed  to  fig. 
(1)  and  eqn.  (28)  in  the  paper  by  Rehm  and  Baum  [3].  Two  similar  examples  can 
be  given,  the  first  being  the  forced  flow  in  the  initial  stages  of  a fire. 

For  this  case,  expansion  in  the  room  of  fire  origin  forces  air  into  an 
adjacent  compartment.  For  this  case,  eqn.  (2)  becomes  (P'  ~ ~ P3) 


1_  ^ 

P'  dt 


1 

(3  - 1)V 


C C,6T  p 
P d 


3/2 


with  P^  = P^(t=o).  If  the  asymptotic  approximation  were  correct  then  the 
right  hand  side  would  be  constant  for  an  e-foldlng  time  for  P'.  Integrating 
this  equation,  we  find  a characteristic  time  for  the  pressure  equation  is 

V 


where  V is  the  volume  of  the  compartment  of  interest,  A is  the  vent  opening 
and  6T  is  the  temperature  difference  between  the  two  compartments.  For 
typical  initial  conditions  we  obtain  a time  constant  about  1 second,  but  it 
can  vary  over  a wide  range.  This  is  certainly  not  less  than  the  typical  time 
step  used  for  the  algebraic  solver  in  fire  models.  The  other  extreme  would  be 
a near  steady  state  regime  of  equal  inflow  and  outflow.  In  this  case,  the 
equation  reduces  to 


dt 


= ,0^  A |.l/2  (6T)^'^^ 

(RT)^ 


,3/2 
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where  R is  the  gas  constant.  For  typical  values  of  the  parameters,  this 
equation  has  a time  constant  of  about  0.10  seconds,  once  again  being  within 
the  range  of  times  used  in  solving  the  other  conservation  equations. 


5.  CONCLUSIONS 


We  have  tried  to  indicate  where  the  problems  arise  in  solving  the 
conservation  equation  used  in  predicting  fire  growth  and  smoke  spread.  The 
type  of  predictive  equations  used  and  the  primary  driving  term  have  been 
discussed  in  detail.  Finally,  a sample  calculation  has  been  presented  to  show 
where  the  problems  actually  exist  for  equation  solvers  and  the  reason  that  an 
ODE  solver,  in  general,  has  an  easier  job  in  threading  its  way  through  the 
phase  space  thicket  than  does  an  algebraic  solver.  Experience  has  shown  [4] 
that  this  method  yields  a reduction  in  computing  time  of  a factor  of  two  to  an 
order  of  magnitude  over  the  time  required  for  a model  wherein  the  asymptotic 
approximation  has  been  made. 
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